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ABSTRACT 


We  combine  Gaussian  wave  packets  and  the  coupled  channel 
method  to  develop  a  theory  of  diffraction  and  rotational 
excitation  by  collision  with  surfaces.  This  improves  our  previous 
work  on’  diffraction  since  it  eliminates  the  mean  trajectory 
approximation:  it  also  extends  Heller's  work  to  problems  in  which 
the  dynamics  require  the  creation  of  new  packets  which  must  be 
coupled  to  each  other  as  they  are  propagated  through  the 
interaction  region.  The  approximations  involved  in  the  above 
Gaussian  wave  packet  approach  can  be  removed  by  using  extending  a 
method  proposed  by  Fleck,  Morris  and  Feit,  which  propagates  the 
Gaussian  wave  function  exactly  and  efficiently. 


I. 


INTRODUCTION 


The  analysis  of  various  H.  surface  scattering  experl - 
1-15  ^ 

ments  requires  the  use  of  quantum  mechanics  In  describing  the 
rotational  motion  and  the  translation  of  the  center  of  mass. 
"Exact"  coupled  channel  calculations  are  possible  only  for  low 

1 1  fi 

Incident  kinetic  energy.  However,  even  when  feasible,  such 

calculations  are  tedious  and  perhaps  Insufficiently  descriptive  of 

the  underlying  physical  processes.  As  a  result  much  work  has  been 

done  to  develop  simpler  and  hopefully  more  Illuminating 

19—34 

approximate  procedures . 

In  this  paper  we  present  an  extension  of  our  previous 
35 

work  In  which  we  used  a  Gaussian  wave  packet  (GWP)  mean 
trajectory  approximation  (MTA)  method  to  calculate  the  diffraction 
and  the  rotational  excitation  of  colliding  with  a  solid 
surface.  The  GWP-MTA  theory  uses  a  wave  function  of  the  form 


♦  (1?,e,*;t) 


N 

=  2 
a»l 


G^(^;t) 


n 

2  c  (t)  Y,  (s,(|»)exp{-le.t/h} , 
1=1  ^  ^ 


(I.l) 


where  ^  Is  the  center  of  mass  position  and  e  and  ^  are  the  polar 
and  azymuthal  angles  describing  the  orientation  of  the  molecular 
axes  In  the  coordinate  system  shown  In  Fig.  1.  The  functions 
Y^(e,^)  are  spherical  harmonics  labelled  by  the  binary  Index  1  = 
(l,m^),  and  G^($,t)  are  Gaussian  wave  packets.  The  experimental 
conditions  are  such  that  the  vibrational  energy  of  H2  exceeds  all 
other  energies  In  the  system,  so  we  can  consider  to  be  a  rigid 
rotor . 

Since  we  are  Interested  In  diffraction  we  must  ensure  that 

the  initial  state  of  the  center  of  mass  Is  a  plane  wave.  This  is 
35—39 

achieved  by  placing  the  packets  G^,  a=l ,  . . . ,M  on  a  M  point 

square  grid  which  covers  the  surface  unit  cell,  and  by  choosing 
the  parameters  in  each  G  so  that  for  ^  within  the  unit  cell,  S 
G^  (^)  coincides  with  the  incident  plane  wave.  The  translational 
symmetry  of  the  surface  allows  us  to  construct  the  result  of 


2 


scattering  by  the  whole  surface,  from  the  results  of  scattering 

35—39 

the  packets  by  an  unit  cell. 

To  explain  some  of  the  improvements  contained  in  GWP-MTA  we 

40—50 

compare  it  with  the  customary  classical  MTA  (denoted  CMTA) 

which  has  been  applied  fairly  successfully  to  surface  science 
28  29  45 

problems.  '  '  When  applied  to  scattering  CMTA  replaces  in 

the  Hamiltonian  the  quantum  variable  ^  with  the  "classical" 

trajectory  5(t),  and  uses  the  wave  function  4»(e,^)  »  Z  c^(t) 

Y^(e,<^)exp{-le^t/h} .  The  classical  trajectory  ^(t)  is  obtained  by 

solving  Newton's  equation  with  the  mean  potential  S  Z  c^(t)Cj(t) 

<Y^ 1 V(^{ t) ) 1 Yj>  where  V(^(t))  Is  the  surface-molecule  Interaction 

energy.  By  using  a  classical  trajectory  MTA  eliminates  all 

quantvim  effects  from  the  motion  of  the  center  of  mass,  except  for 

those  contained  in  the  computation  of  the  mean  force.  The  GWP-MTA 

theory  propagates  a  Gaussian  wave  packet  on  the  same  potential 

energy  as  CMTA.^®'®®  The  use  of  a  wave  packet  provides  a  fully 

quantum  mechanical  description  (albeit  a  simplified  one)  of  the 

center  of  mass  motion.  The  resulting  theory  describes  well  such  a 

dramatic  quantum  effect  as  diffraction  and  also  gives  a  good 

35 

description  of  the  rotational  excitation  probabilities. 
Furthermore,  GWP-MTA  is  computationally  cheaper  than  its  classical 
counterpart  since  the  number  of  GWP's  needed  to  describe  the 
quantum  scattering  process  is  substantially  smaller  than  the 
number  of  classical  trajectories  required  by  CMTA.  This  can 
generate  substantial  savings  in  computer  time  since  the  expense 
per  GWP  is  roughly  twice  the  expense  needed  to  propagate  a 
classical  trajectory.  This  can  be  understood  metaphorically,  by 
thinking  of  each  packet  as  a  bundle  of  classical  trajectories, 
which  are  generated  at  once  by  propagating  one  packet. 

Unfortunately  GWP-MTA  shares  with  CMTA  a  shortcoming  whose 
removal  is  the  subject  of  the  present  paper.  To  understand  both 
the  origin  of  this  shortcoming  as  well  as  its  removal  by  the 
method  proposed  here  it  is  useful  to  compare  the  GWP-MTA  wave 
function  (II.l)  with  its  proposed  replacement,  which  is 


(1.2) 
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♦  t) 


M  n 
*  Z  Z 
a*l  1*1 


Y^(e.^)exp{-le^t/h} . 


The  GWP-MTA  wave  function  describes  the  center  of  mass 
motion  with  one  Gaussian  per  grid  point.  Superficially  this  is  in 
agreement  with  out  intuitive  belief  that  a  paclcet  represents  a 
localized  "corpuscule"  (as  opposed  to  a  wave)  and  therefore  we 
must  have  one  trajectory  per  corpuscule.  It  is  however  incorrect 
to  apply  this  "Newtonian"  notion  to  a  system  that  has  discrete 
internal  states  (i.e.  rotations)  which  can  be  excited  during 
collision.  This  makes  the  classical  motion  of  the  center  of  mass 
rather  non-Newtonian,  since  the  quantum  excitation  of  the  internal 
motion  affects  the  center  of  mass  motion  (at  least  through  the 
conservation  laws).  Thus  the  excitation  of  an  internal  state 
requires  the  appearance  of  a  new  center  of  mass  trajectory  whose 
translational  energy  is  equal  to  the  incident  one  minus  the  energy 
of  the  internal  excitation.  Therefore,  a  correct  description  of 
center  of  mass  motion  requires  one  trajectory  for  each  final 
rotational  state.  Within  the  GWP  approach  this  can  be  achieved  by 
using  the  wave  function  (1.2)  which  has  (at  each  grid  point)  one 
packet  for  each  rotational  state. 

It  is  now  useful  to  contrast  the  behavior  of  the  packets  in 
these  two  theories.  The  unique  packet  used  in  GWP-MTA  moves  on  a 
mean,  rotatlonally  averaged  potential  surface.  A  time  of  flight 
(TOP)  measurement  applied  to  this  theory  gives  one  peak  in  the 
aomentiuB  distribution  whose  position  is  determined  by  the  fact 
that  the  kinetic  energy  lost  by  the  center  of  mass  motion  equals 
the  total  energy  taken  up  by  all  the  rotations;  the  same  TOP 
measurement  on  a  system  described  by  the  multiple  GWP  theory 
(MGWP)  leads  to  one  peak  for  each  rotational  state.  In  MGWP  the 
post-collision  rotational  distribution  is  Imprinted  in  the  TOP 
spectra  and  with  sufficient  resolution  one  could  measure  the 
rotational  energies  by  doing  TOP  measurements.  Another  way  of 
pointing  out  differences  between  the  two  methods  is  to  examine  the 
results  predicted  for  rotatlonally  selected  final  state 


measurements.  GWP-MTA  gives  the  same  diffraction  spectrum  for 
each  rotational  state,  while  MGWP  gives  different  diffraction 
patterns  for  each  final  rotational  state.  Furthermore  for  a 
rotationally  selected,  angle  resolved,  TOP  measurement  GWP-MTA 
gives  the  same  results  for  each  state  Y^,  while  MGWP  predicts 
different  results  for  different  Y^ . 

We  emphasize  that  while  the  conceptual  and  qualitative 

Improvements  brought  about  by  MGWP  method  are  Interesting  and 

necessary,  one  should  not  view  the  use  of  GWP-MTA  with  exagerated 

alarm.  The  measurements  required  to  discriminate  between  the  two 

methods  are  possible  but  very  tedious.  Less  ambitious  ' (but  still 

difficult)  experimental  work,  such  as  a  measurement  of  the 

rotational  distribution  with  a  modest  angular  resolution  and 

without  TOP,  is  likely  to  be  well  described  by  GWP-MTA,  since  the 

method  of  measurement  performs  experimentally  the  kind  of 

averaging  that  GWP-MTA  does  theoretically.  Diffraction 

measurements  with  modest  angular  resolution  and  no  anlysis  of  the 

rotational  state  are  also  well  described  by  GWP-MTA  except  for 

those  situations  when  the  angular  resolution  is  sufficient  to 

resolve  the  diffraction  peaks  due  to  molecules  that  are 

2  6 

rotationally  excited.  '  Such  peaks  are  averaged  by  the  MTA 
theory  together  with  the  rotationally  elastic  parent  peak  and  will 
be  absent  in  the  predictions  of  the  theory. 

Within  the  existing  GWP  methodology®^”^®  the  development  of 
the  MGWP  scheme  requires  the  solution  of  several  technical 
problems.  The  customary  propagation  scheme  assumes  that  the 
number  of  packets  is  conserved  in  time,  while  in  the  present 
theory  we  must  start  with  one  packet  and  emerge  with  as  many 
packets  as  many  open  rotational  channels.  Thus  we  must  find  a 
method  for  generating  new  packets  as  the  incident  packet  enters 
the  collision  zone. 

Furthermore,  the  existing  applications  of  the  GWP  method  to 
3  6“  3  9 

diffraction  has  used  what  we  call  the  simplest  Heller  method 

(SHM),  which  assumes  that  the  packets  can  be  propagated 
independently.  This  assumption  cannot  be  made  in  the  present 


5 


problem  for  the  Gausslans  .  1=1 »  m  because  the  rotational 

populations  are  established  exclusively  through  the  coupling 
between  the  packets  having  the  same  a  but  different  1-s.  Neglect 
of  this  coupling  (as  in  SHM)  would  suppress  rotational  excitation 
from  the  theory. 

The  main  contribution  of  the  present  paper  to  the  GWP 
technology  is  the  development  of  a  propagation  scheme  in  which 
Gausslans  are  created  in  the  process  of  propagation  as  needed,  and 
evolve  by  interacting  with  each  other. 

The  difficulties  encountered  by  a  classical  or  semi- 
classical  propagation  scheme  applied  to  a  subset  of  degrees  of 
freedom  which  are  coupled  to  a  strongly  quantum  subset,  are  well 

71-79 

known  and  much  labor  and  ingenuity  has  been  devoted  to  their 

solution.  The  present  MGWP  approach  is  a  new  procedure  to  attack 
this  old  problem.  Since  the  space  does  not  permit  a  detailed 

71-79 

comparison  between  the  present  and  the  earlier  work,  we 

confine  ourselves  to  listing  those  features  that  make  us  hope  that 

the  method  developed  here  will  be  useful.  First,  the  MGWP  does  not 

require  root  searching,  classical  calculations  with  double  ended 

boundary  conditions,  or  self-consistent  solutions  of  integro- 

dlfferentlal  trajectory  equations,  it  is  not  confined  to  one 

dimension,  and  has  no  difficulties  with  the  turning  points. 

Second,  the  GWP  equations  of  motion  are  almost  as  simple  and 

sometimes  less  laborious  than  the  classical  ones.  The  GWP 

scheme  blends  easily  and  naturally  with  classical  mechanics  so 

that  quantum  scattering  calculations  from  a  classically  moving 

39 

lattice  are  possible.  Finally  the  method  lends  Itself  to  simple 

classical  like.  Intuitive  interpretation  of  dynamics. 

In  the  context  of  surface  science  GWP  methods  are  relatively 

new;  the  existing  calculations  show  that  they  are  reasonably 
35—39 

accurate.  In  diffraction  problems  they  can  be  applied  at 

kinetic  energies  at  which  coupled  channel  calculations  are 

prohibitive.  They  can  be  easily  used  to  calculate  scattering  by 

37 

disordered  surfaces,  a  problem  for  which  the  traditional  quantum 
methods  would  have  extreme  difficulties. 


The  single  most  important  drawback  of  the  GWP  method  is  that 

it  is  an  "ill-defined"  approximation,  in  the  sense  that  it  lacks  a 

precise  validity  condition,  or  a  convergence  scheme  which  insures 

the  achievement  of  greater  accuracy  with  increased  labor.  This 

cannot  be  done  by  increasing  the  number  of  packets;  in  fact  it 

sometimes  happens  that  an  increase  in  the  number  of  Gaussians 

67  69 

leads  to  overcompleteness  and  worsens  the  accuracy. 

Practical  experience  indicates  that  one  should  expect  good  results 

with  little  effort  for  short  time  dynamical  problems  involving 

spatially  localized  quantum  degrees  of  freedom,  which  are  nearly 

seml-classlcal .  Recent  work  by  Sawada,  Heather,  Jackson  and 
67 

Metlu  shows  that  for  some  problems  accurate  long  time  results 
can  also  be  obtained  at  the  expense  of  the  simplicity  of  the 
propagation  scheme. 


II.  THE  MULTIPLE  GAUSSIANS  THEORY  OP  H2  SCATTERING 
II. 1.  The  Hamiltonian 

He  use  the  Hamiltonian 


H  »  -{h^/2M)y^ 
R 


+  {2u|rr) 


^2 


(II. 1) 


where  ^  is  the  center  of  mass  position,  r^  and  r^  are  the 
coordinates  of  the  two  atoms  forming  the  diatomic,  M=2M  and  |j=M/2 
are  the  total  and  the  reduced  mass  of  the  diatomic,  and  t,  is  the 
angular  momentum  operator.  The  coordinate  system  is  shown  in  Fig 
1.  The  interaction  energy  with  the  surface  is 


V(ri,r2)  -S  D  { exp[ -2a^  ( z  ^-z^)  ]-2exp[ -ot^  ( z^-z^ )  ]  } 

(II/2) 

-PjD^exp[-2aj^(z^-z^)  ]  [cos(2-rrXj|^/Cj^)  +  cos  (  2TTy^/c2 )  ]  . 

where  c^^  and  C2  are  the  lattice  constants,  and  rj^=(x^,y^,z^)  . 

This  represents  the  interaction  with  a  corrugated  surface  and  the 
forces  acting  on  the  molecule  depend  on  both  the  polar  and  the 
azymuthal  angle  describing  the  orientation  of  the  molecular  axis 
(Fig,  1).  To  obtain  the  equations  of  motion  for  the  nuclear  wave 
fiinctlons  G^^  we  Insert  the  wave  function  (1.2)  in  the  time 
dependent  Schrodlnger  equation,  use 

(2ulr|2)”l  L2Y^(e,4,)  =  e^Y^(e,(|>)  ,  (11.3) 


multiply  with  Yj(e,^)  sineded<>,  and  integrate  over  angles.  This 
leads  to 


Cih 


at 


•(h2/2M)72]Z  G  (^;t) 
^  a  aj 


=  5  Vj^(^)exp(-i(ej-e^)t/h) (ZG^. ) . 


with 

v,,(S) 


TT  211  * 

;  sined©  s  d(|)  Y  .{©.((»)V(r,  ,r,)Y.  (©,4))  . 


(II. 5) 


II . 2  The  Equations  of  Motion 


II. 2. A.  General  Remarks 

To  use  the  GWP  method  for  solving  the  equations  (II. 4)  we 
assume  that  G^^  are  wave  packet  of  the  form 


G^^(^;t)  =  exp[(i/h){(^-^^^(t))-*t^.(t).(^-^„i(t)) 


(II. 6) 


Here  and  are  the  expectation  values  of  the  position 

and  momentum  operators  in  the  state  G  A  .(t)  is  a  3x3  comolex 

ai  ai 

matrix  which  gives  the  width  of  the  packet  and  a  space  dependent 
phase  and  is  a  complex  function  of  time  contributing  to  the 

phase  and  the  amplitude  of  the  state 

The  central  idea  in  Heller's  work  is  that  the  time  evolution 
of  the  state  G^^  is  given  by  the  evolution  of  the  parameters 
t),  In  the  simplest 

implementation  of  this  idea  (which  we  call  the  simplest  Heller 
method  (SHM))  it  is  assumed  that:  (1)  the  Gaussians  are  narrow 
throughout  the  collision  so  that  we  can  expand  thte  potential  in 
the  Schrodinger  equation  (II. 4)  in  power  series  around  the 
instantaneous  position  ^^^(t)  of  the  packet  and  retain  che 
quadradic  part  only  (the  local  harmonic  approximation  (LHA));  and 
that  (2)  each  Gaussian  can  be  propagated  independently  (the 
Independent  Gaussian  approximation  (IGA)).  The  shortcomings  of 

6  7-  6  .9 

SHM  were  pointed  out  in  our  work  ate  well  as  the  work  of  Skodje  and 
Truhlar®®  and  Thlrumalal,  Bruskin  and  Berne.®®  Ways  of  improving 
SHM,  by  removing  the  two  approximations  presented  above  were 
proposed  by  Sawada,  Heather,  Jackson  and  Metiu®^  and  by 
Heller.®^'®® 

In  this  paper  we  use  the  minimum  error  method  (MEM)  which 
couples  the  Gaussians  and  makes  no  assumption  concerning  their 
width.  The  only  approximation  is  that  G  .(l^;t)  maintain  their 
Gaussian  form  throughout  the  collision. 

Our  numerical  experience  with  the  coupled  Gaussian  MEM 
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theory  Indicates  that  the  coupling  terms  vary  on  a  faster  time 

scale  than  the  other  terms  appearing  In  the  MEM  equations.  The 

presence  of  these  terms  slows  down  considerably  the  integration 

program  and  therefore  it  is  desirable  to  neglect  them  whenever  the 

penalty  is  not  too  severe.  In  the  diffraction  calculations  carried 

36  37 

out  by  Drohlshagen  and  Heller  '  (for  atoms)  and  by  Jackson  and 
38  3d  35 

Metiu  (for  atoms  '  and  diatomlcs  )  it  was  demonstrated  that 
satisfactory  results  can  be  obtained  if  the  Gausslans  used  to 
construct  the  Incident  plane  wave  are  propagated  Independently. 

In  the  present  context  this  suggests  that  we  can  assume  that 
Gpj  etc  are  decoupled  when  We  cannot  however  neglect  the 

coupling  between  the  packets  G^^  1=1,  ....  n,  since  this  coupling 
is  the  driving  force  for  the  rotational  excitation  of  the 
particle. 

II. 2. B.  The  Equations  of  Motion 

To  generate  the  MEM  equations  of  motion  for  the  Gaussian 
67 

parameters  we  define  the  error  made  by  using  the  Equation  (1.2) 


as : 


E  =  J 


(h^/2M)7^  +  lh|^)  G^j(^;t)-  Z  ^  (i5)exp{  1  (  e^-e^ )  t/h) 
R  i  ®  1 


,2 

G^^(^;t)  d^ 


Since  the  error  made  at  time  t  is  a  function  of  A  ,(t),  R  ,(t), 
j,  ai  ai 

Po^l(t)  and  we  determine  these  quantities  so  as  to  minimize 

the  error.  This  generates  first  order,  non-linear,  ordinary 
differential  equations®^  for 


^ai 


(t)  : 


■  S  (II. 7a) 


(II. 7b) 


(11.70 


“J  “J  ”  ^^ooo^aj  ^^ooo’aj 


F^"’=  Z  Xd]5(^-^  (t)"  G*j(^;t)  G^^(^;t)V  (5)e'^^®i'^j 


(II. 7d) 
,-eOt/h 


(II. 8) 

and 

*TrJ^’=  Xd^(^-^^j(t)  .  (S-^^j(t)  |G^^.(i^;t)  l^  .  (II. 9) 

^(O 

Note  that  F^j  is  a  scalar,  F^^  is  a  vector  with  the  components 

(F^J^x  -  I  /d5(X-X^j(t)G*jG^^Vj^exp(i(ej-e^)t/h)  ,etc. 

2 ) 

and  F^j  is  a  dyadic  with  the  components 
'^i3’’xy  =  JM*(J'-x„j(«)>(V-Y^j(t)) 

^exp(-i(ej-e^)t/h) ,  etc.. 

The  complex  matrix  S  .  and  the  vectors  I  .  and  X  .  are 

Qtj  «j  0(j 

defined  in  the  Appendix,  where  we  also  give  an  outline  of  the 
derivation  of  the  equations  presented  above. 
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II . 3  The  Choice  of  Initial  Condition 

In  order  to  solve  the  equations  (II. 7)  we  must  choose 
initial  values  for  the  parameters  A^^{0), 

As  discussed  often  in  literature  the  GWP  method  does  not  provide 

an  unique  and  fully  justifiable  method  for  specifying  these 

values.  However  in  the  particular  case  of  atom  or  diatom 

diffraction,  prescriptions  that  work  well  have  been  proposed  and 

35^39 

we  see  no  reason  to  modify  them.  Thus  if  the  initial 

rotational  state  is  we  take  for  j?«i  and  define  the 

parameters  of  such  that  is  s  plane  wave  when^  takes 

values  within  the  unit  cell  of  the  solid  surface.  The  way  to  do 

35—39 

this  was  discussed  in  several  papers. 

While  this  provides  initial  conditions  for  the  equations 
propagating  the  parameters  in  ,  it  says  nothing  about  the 
initial  conditions  for  There  is  a  physical  reason  for 

this:  the  initial  state  Z  G^^  is  prepared  by  the  experimentalist, 

while  Z  G^j  is  generated  by  the  collision  process.  The  available 
GWP  procedures  have  no  provisions  for  the  birth  of  new  packets . 

To  circumvent  this  difficulty  we  propose  the  procedure  described 
below. 

We  consider  the  case  in  which  the  incident  molecules  are  in 
state  and  the  incident  state  of  the  center  of  mass  is 
constructed  with  the  packets  ot=l ,  ....  M.  Let  us  denote  by 

t^  the  last  time  when  we  can  still  assume  that  the  packets  G 
l?!k  have  zero  amplitude.  At  that  time  the  packets  satisfy  the 
equations 


(Ih  +  (h^/2M)V^}  G^j^(5;t)  =  Vj^j^(^)G^j^(R;t) 


( II . 10a) 


(ih  +(h^/2M)7|)G^^(5;t)  =  V^j^(^)exp{-l(ej^-e^)t/h}G^j^(^;t) 


l?Jk,  1=1,2, 


II . 10b) 


According  to  these  equations  a  new  packet  G  ^(^;t)  is  generated 
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when  the  initial  packet  starts  overlapping  with 

To  create  the  new  packet  solve  Eq(II.lOb)  for 

a  short  time  t.  This  gives 

«  -(iT/h)  {M/ih2TTT)^'^  J*  d$' 

— W 

exp{  (  i/h)  [  (M/2T)  )  2]  )Gj^(^'  ;  t^) 

exp{-{i/h)  (ej^-e^)t^}  .  (II. 11) 

A  rather  tedious  calculation  shows  that  for  the  surface- 

molecule  potential  (II. 2)  the  quantity  )  »  defined  by  Eq. 

(II.S),  is  a  Gaussian.  Therefore  the  integrand  in  Eq.  (II.  11)  is 

a  product  of  Gaussians  and  therefore  the  integration  gives  a 

Gaussian  for  If  ^Ik^^'^  ®  complicated  non-Gaussian  form 

the  Integral  can  be  very  accurately  performed  (since  t  is 

arbitrarily  small)  by  using  the  stationary  phase  approximation  and 

this  also  gives  a  Gaussian  result. 

It  should  be  noted  that  the  above  procedure  is  nothing  else 

but  perturbation  theory  with  respect  to  followed  by  an 

asymptotic  evaluation  of  an  Integral  by  using  the  fact  that  the 

1/2 

inverse  length  (M/2Th)  is  very  large  compared  to  all  other 
length  scales  in  the  problem  (i.e.  the  width  of  G^j^  and  the  length 
scale  over  which  Vj^j^(^)  changes).  Because  we  control  the 
magnitude  of  t  this  procedure  is  essentially  exact. 

80 

Our  numerical  experience  with  charge  transfer  problems  shows 

that  it  is  possible  to  generate  all  packets  at  once,  even  if  for 

some  of  them  the  term  is  very  small  due  to  poor  overlap. 

The  outcome  is  that  the  new  packet  G  has  a  much  smaller  initial 

ocm 

amplitude  than  the  other  new  packets,  as  it  should.  Subsequent 
propagation  will  Increase  that  amplitude  as  Increases. 

A  serious  danger  in  using  such  a  method  is  that  the  results 
may  depend  on  the  moment  t^  chosen  for  new  packet  generation. 

There  is  some  numerical  experience  regarding  this  possibility, 


k& 


which  indicates  that  the  results  are  remarkably  stable  with 
respect  to  this  choice. 

II. 4  A  Simplified  Model  and  A  Qualitative  Discussion  of  the 


Theory 

The  notation  used  so  far  gives  a  compact  representation  of 
the  equations  of  motion  (II. 7)  for  the  Gaussian  parameters,  but 
obscures  to  some  extent  their  physical  meaning.  In  this  Section 
we  discuss  a  limiting  case  which  is  designed  to  simplify  the 
cpmputatlonal  scheme  and  bring  out  its  physical  content.  The 
approximations  made  are  sufficiently  reasonable  to  give  us  hope 
that  the  simplified  scheme  might  also  be  accurate  enough  to  be 
computationally  useful. 

As  Heller  pointed  out,  in  the  simplest  version  of  his  method 
(SHM)  ^(t)  and  ?(t))  follow  classical  equations  of  motion  and  the 
phase  Re7(t)  is  essentially  the  classical  action.  These  results 
are  obtained  only  if  the  Gaussians  are  narrow  (LHA)  and  decoupled 
(IGA)  . 
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The  MEM  version  of  the  theory  does  not  make  these 
simplifying  assumptions  and  as  a  result  loses  the  simple  features 
mentioned  above.  Nevertheless,  it  is  still  useful  to  think  of  the 
equation  of  motion  for  and  as  evolving  according  to  a 
mechanics  which  is  similar  to,  but  more  complicated  than  the 
classical  one;  for  lack  of  a  better  term  we  call  this  a 
corpuscular  mechanics  (we  have  also  used  in  our  work  the  term 
pseudo-classical).  When  the  MEM  theory  uses  only  one  packet  G^^ 
the  force  acting  on  the  center  of  the  packet  is  not  -3V.(R  ,)/3R  . 

^  ,  1  ai  Qti 

(as  in  the  simple  Heller  method)  but  -3/3R  .  /  G  ,  V.G  ,dR.  Since 

n  oci  al  1  oci  ^ 

the  quantity  |G^^|  appearing  in  this  integral  depends  on  ImA(t), 

the  above  force  depends  on  time  and  the  motion  of  the  center  of 

the  packet  is  not  conservative.  Various  interesting  features  of 

this  potential  and  the  resulting  "corpuscular  mechanics"  were 
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discussed  and  exemplified  by  Heather  and  Metiu.  The  use  of 
several  coupled  Gaussians  further  complicates  the  picture  since 
the  motion  of  their  centers  is  now  coupled  through  terms  that  have 
no  classical  analog  and  which  are  turned  on  by  the  overlap  between 


the  packets  (i.e.  they  depend  on  the  quantity  J*  dR(R-Rj.) 

G;6jV,j{R)). 

In  what  follows  we  analyze  the  meaning  of  the  MEM  equations 
(II. 7)  by  making  some  of  the  approximations  used  by  the  SHM.^^ 
First  we  emphasize  again  that  while  in  previous  applications  of 
MEM  we  coupled  the  Gaussians  in  order  to  increase  the  accuracy  of 
the  calculations,  in  the  present  work  the  coupling  between  the 
Gaussians  ,  i=l,...,n  with  the  same  a  is  essential  on  physical 
grounds.  It  is  this  coupling  that  allows  the  amplitudes  of 
various  Gaussians  to  vary  and  give  us  the  final  rotational 
distribution.  The  probability  that  a  particle  is  in  the 
rotational  state  is  proportional  to  g  *^al^ai  depends 

on  Imy^^(t)  and  det  ( ImA^^  ( t ) )  .  Thusy^compute  this  probability 
correctly  it  is  essential  that  the  coupling  between  the  Gaussians 
G^j^ ,  1*1,  ...,  n  is  present  in  the  equation  of  motion  for  YQjj(t) 
and  however,  we  might  expect  that  the  coupling  may  be 

neglected, ' without  causing  a  dramatic  qualitative  deterioration  of 
the  results,  in  the  equations  of  motion  for  and  • 

If  we  begin  with  the  Equation  (II. 7a)  for  ^  .  we  find  that 
the  coupling  between  packets  enters  through  Im  .  This  quantity 
can  be  written  as  (see  Eq.  (II.9)). 


z  Im?,  . 
^  J  1  ^  ^ 


with 


(11.12) 


?^j«Xd^(^-5^j(t) )G*j(R;t)Vj^(^)G^^(^;t)exp(-i(e^-e^)t/M 


(11.13) 

The  diagonal  term  ?jj(t)  is  real  and  therefore  does  not  contribute 
to  Im  j  •  We  are  thus  left  with  the  off-diagonal  terms  . 
only.  These  represent  the  effect  of  the  coupling  of  G  .to  the 

J 

other  Gaussians  G^^ ,  i;^ j .  Since  we  argued  that  in  the  equation 
for  this  coupling  is^^^ot  physically  essential  we  neglect  it, 
and  therefore  can  take  ?  .  0 . 


15 


This  assumption  is  also  reasonable  in  view  of  our  past 
experience  with  coupled  Gaussian  calculations.  We  have  found  that 


the  time  evolution  of  the  terms  G  .G  4  is  much  more  rapid  than 

that  of  or  ^  because  of  rapid  changes  in  the 

difference  between  the  phases  of  G  .  and  G  . .  Thus  the  terms 

^  ai  otj 

^ijji^j,  oscillate  very  rapidly  and  their  effect  on  the  evolution 


of  may  average  to  zero  when  the  equation  of  motion  for  is 

integrated;  furthermore  the  phase  exp{ i { e .-e  . ) t/h}  can  play  the 

j  « 

.  .  is  non-zero  only  when  G  . ,  V, ^  and  G. 


- -  - - -  -  - - 

overlap  which  means  that  is  practically  zero  at  least  at  some 

times  during  collision. 


If  these  arguments  are  accepted  we  can  replace  (11.7a)  with 

(11.14) 


^  .  a:  ?  ./M 
aj  cxy 


which  implies  that  the  center  of  each  packet  G^j  moves  like  a 
classical  particle  with  moment\im  * 

A  similar  discussion  can  be  made  to  simplify  Eq.(II.7b) 
(giving  the  evolution  of  )  by  neglecting  when  ixj .  The  use 
of  Eq.  (IZ.14)  in  (II. 7b)  removes  the  last  term  in  the  latter 
equation.  If  we  further  assume  that  the  Gaussian  G^^  is  narrow 
throughout  collision  we  can  expand  V^j(^)  appearing  in  the 
integrand  of  (see  Eq.  (11.13))  and  obtain 


it 


Sd^  (^-I5^j(t)  >(?-i^^j(t)  )  lG^^(R;t) 


The  last  term  follows  from  Eq.  (II. 9).  Using  these  approximations 
and  Eq.  (11.12)  in  Eq.  (II. 7. b)  leads  to 


(t) 


-3V(^^j(t)  )/3l5^j(t)  . 


( II . 15a) 


If  we  don't  assume  a  narrow  Gaussian  we  have 
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{11.15b) 

Thus  the  use  of  the  Independent  Gausslans  approximation  (IGA)  in 

(II. 7a)  and  (II. 7b)  leads  to  (11.14)  and  (II. 15b);  the  additional 

assumption  that  the  potential  energy  is  locally  harmonic  (LHA) 

simplifies  (Il.lSb)  and  gives  (II. 15a).  If  both  LHA  and  IGA  are 

used  each  packet  center  moves  classically  on  the  potential 

Vjj(^^j(t)),  which  is  the  expectation  value  of  the  molecule 

surface  potential  energy  when  the  molecule  is  in  the  rotational 

state  .  If  LHA  is  abandoned  the  motion ^^esembles  superficially 

the  classical  one.  However  the  force  -[M  is  time 

«J  jj 

dependent  and  the  motion  of  the  packet  is  non-conservative.  As 

6d 

shown  by  Heather  and  Metiu  '  this  is  a  necessary  feature  which 
ensures  the  conservation  of  energy;  if  LHA  and  IGA  are  used  the 
"classical  energy"  +  V(lS^^(t))  is  conserved,  but  the 

quantum  one  isn't. 

We  emphasize  that  the  use  of  IGA  to  decouple  the  motion  of 
the  centers  of  the  packets  is  not  likely  to  lead  to  qualitative 
errors  for  short  collision  times,  and  simplifies  considerably  the 
propagation  scheme  through  the  elimination  of  the  rapidly  varying 
coupling  terms.  There  is,  as  yet,  no  proof  that  the  errors  made 
by  using  these  simplified  equations  are  small  and  the  usefulness 
of  these  equations  remains  to  be  tested.  In  support  of  this 
simplified  theory  we  note  that  it  is  superior  to  all  classical 
trajectory  methods  which  use  one  trajectory  only  since  it  has,  as 
is  conceptually  required,  one  trajectory  for  each  rotational 
state.  On  conceptual  grounds  the  theory  also  compares  favorably 
with  the  well  known  Preston-Tully  (PT)  method,^®  while  being  less 
demanding  numerically.  Like  PT  we  have  multiple  trajectories,  on« 
for  each  rotational ly  averaged  energy  surface.  Unlike  PT  the  GWP 
procedure  -  even  in  its  simplest  form  discussed  in  this  Section- 
describes  the  center  of  mass  motion  by  using  a  nuclear  wave 
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function.  The  trajectories  are  only  a  simplifying  device  for 

propagating  these  wave  functions.  As  a  result  the  present  theory 

uses  probability  amplitudes  while  PT  attaches  to  each  state  a 

probability,  ignoring  the  superposition  principle  and  losing 

interference  effect.  Thus,  the  present  theory  describes 

interference  dominated  phenomena  -  such  as  diffraction  -  rather 

well  while  PT  cannot  describe  them  at  all.  Furthermore  the  use  of 

GWP’s  (with  (11.15b))  maintains  Heisenberg  principle  in  the  theory 

and  this  should  improve  the  dynamics,  especially  if  zero  point 

motion  is  important.  Finally  we  compute  observables  describing 

center  of  mass  motion  by  using  the  rules  of  quantum  mechanics 

(i.e.  wave  functions,  operators,  matrix  elements,  etc.)  while  the 

PT  method  is  confined  to  classical  rules.  Thus,  for  example,  all 

*  2 

effects  of  quantum  fluctuations  (e.g.  the  fact  that  <G|P  \G  x 
*  2 

<G|P|G>  )  are  lost  in  the  PT  method.  Nevertheless  while 
conceptual  improvements  are  pleasing,  a  direct  numerical 
comparison  between  MGWP  and  PT  is  required  to  test  whether  such 
improvements  have  any  practically  useful  consequences. 
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III.  DISCUSSION 

III.l  Remarks  regarding  the  errors  made  by  MGWP 

The  errors  in  MGWP  are  made  because  ( 1 )  we  propagate  the 
packets  originating  from  different  grid  points  independently  and 
(2)  we  force  them  to  maintain  their  Gaussian  form  throughout  the 
collision. 

The  statement  that  it  is  erroneous  to  propagate  i|i  »  I  G  by 
propagating  each  G^  Independently  seems  to  be  at  odds  with  the 
linearity  of  Schrodinger  equation.  Indeed  if  we  apply  to  i|>(0)  =  Z 
G^{0)  the  exact  propagator  U{t)  to  calculate 

♦  (t)  =  U(t)Hi(0)  =,Z  U(t)G^  .  (III.l) 

a 

it  automatically  follows  that 

<t(t) |A|«(t)>  -  Z  <G  (0) |U(t)*AU(t) |G  (0)>  (III. 2) 

1  “  '  “  ♦ 

I  for  any  operator  A.  If  we  take  A°l  and  if  U(t)  » 

0(t)”^,  the  propagation  scheme  will  conserve  the  norm;  if  A“H  and 

U(t)H«HU(t)  the  conservation  of  energy  follows.  These  conclusions 

hold  even  if  U(t)  is  approximate,  as  long  as  it  is  unitary  and 

I  commutes  with  the  H€uiiiltonlan. 

However  the  application  of  the  GWP  method  with  independent 
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Gaussians  to  several  (but  not  all)  examples  shows  that  the  above 

conservation  does  not  hold;  moreover,  if  the  Gaussians  are 

i  coupled,  the  conservation  properties  are  restored. 

This  discrepancy  between  the  conclusions  of  the  familiar 

analysis  presented  above  and  the  numerical  calculation  is  due  to 

the  peculair  nature  of  the  GWP  propagation  method  which  calculates 

I  0(t)G^(0)  by  separately  optimizing  U(t)  for  each  packet  G^.  Thus 

the  result  of  such  a  calculation  is  more  properly  denoted 

U^(t)G^.  As  long  as  it  is  applied  to  G^  the  GWP  propagator  is 

unitary  (l.e.  conserves  the  norm  <G  (0)|G  (0)>)  and  commutes  with 

oc  oc 

I  the  Hamiltonian  (i.e.  <G  (0)|H|G  (0)>  is  conserved).  There  is 


however  no  assurance  that  U  (t)^U  (t)G^  *  nor  that  U  HG„  = 

a  p  p  p  a  P 

HU^Gp.  Because  of  this,  the  quantity 

<i|»(t)  |A|<|>(t)>  =  2  Z  <G  (0)|U  (t)*  OU  (t)|G  (0)>  ,  (III. 3) 
a  p  “  “  P  P 

Is  not  conserved  for  A^l  or  for  A^H.  The  MEM  procedure  does  not 
have  this  shortcoming  since  it  seeks  an  optimum  propagator  for  the 
whole  sum  Z  G  rather  than  for  each  G  independently,  and 

Ct  Ot 

therefore  (II I. 2)  and  Its  consequences  are  valid. 

The  use  of  the  IGA  in  developing  the  MGWP  procedure  Is  based 

on  Its  success  In  previous  diffraction  calculations.  We  have  no 

detailed  understanding  of  the  reasons  for  this  success,  other  than 

the  qualitative  argviment  that  the  IGA-GWP  approximation  gives  each 

packet  a  phase  that  Is  very  similar  to  that  given  by  the  semi- 

8 1 

classical  theory,  which  Is  known  to  work  well  for 
82**85 

diffraction.  Nevertheless  It  is  of  Interest  to  discuss  a 

computational  method  which  removes  the  IGA  approximation  as  well 

as  the  restriction  that  the  packets  maintain  their  Gaussian  form 

during  collision.  This  is  based  on  ideas  advanced  by  Fleck, 
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Morris,  and  Felt  (FMF)  (summarized  in  Section  III. 2)  applied  to 
the  MGWP  theory  of  scattering  (Section  III. 3). 


The  starting  point  of  the  FMF  method  is  similar  to  that  used 
in  path  integral  theory.  The  long  time  propagator  U{T)  is  written 
as 


W(T)  ss  U{T)’^  (III. 4) 

with  nT=T.  The  short  time  propagator  U(t)  has  the  "split"  form 

U(T)=exp{ (iTh^/4M)7|) }exp{-iTV{1?) }exp{ ( iTh^/4M)7|}  (III.  5) 

It  is  easy  to  see  that  U(T)  can  also  be  written  as 

U(T)»exp{  (iTfi^/4M)7|)exp(-iTV(5)  )W(  x  )”"^exp{  iTh^/4M)7|}  , 

(III. 5) 

where 

W(T)  s  exp{ ( iTh^/2M)7|}  exp{-iTV(^)}  (III. 6) 


is  the  more  familiar  expression  for  the  short  time  propagator 
appearing  in  the  discrete  form  of  the  path  integral  formula  for 

3 

U(T) .  The  error  made  by  using  (III. 5)  is  of  order  x  ,  for  each 

time  step  x,  while  the  discrete  path  integral  formula  U(T)  = 

n  2 

W(x)  makes  an  error  of  order  x  at  each  time  step. 

To  explain  the  FMF  algorithm  we  examine  the  computation  of 

the  elementary  step  W(  x  )  1 1((  ( t )  >  .  Using  (III.  6)  and  straightforward 

manipulations  based  on  the  representation  theory  we  can  write 


<1^' 14P(t+x)>  s  <?|;^|W(x)  |.(p(t)>  = 

»  <^^l^^>exp{-(ixh^^J/2M)  }<ic^|5^>exp{-ixV('^^)  } 

■  <^^|^(t)>  .  (III. 7) 


Here  we  have  discretized  the  variables  ^  and  ^  and  used  the  rule 
that  repeated  variables  are  summed  over.  In  the  customary  path 


integral  theory  the  integral  over  is  done  analitically  and 
generates  the  familiar  Gaussian  A  = 

(2TriT/M)  exp{  )  ^/2fiT )  .  We  are  thus  left  with  only  one 

integral,  over  which  must  be  performed  numerically.  This 
amounts  to  a  multiplication  by  the  diagonal  matrix  B^^=6^^exp{- 
iTV(^^)},  followed  by  multiplication  with  the  matrix  .  If  a 
proper  discrete  representation  of  the  function  exp{-iV(^)  }<^|4»> 
requires  a  grid  having  points  (d  is  the  dimension  of  the  vector 
the  calculation  outlined  above  requires  operations  per 

time  step.  This  kind  of  calculation  has  been  performed  by 


Th 1 r tuna laTr ana  Berne  who  used  it  to  evaluate  the  path  integral 

formula  for  the  partition  function,  which  is  the  imaginary  time  ' 

8  8 

verlson  of  our  problem. 

The  PMF  method  does  perform  both  integrals  in  (III. 7)  (i.e. 
the  svims  over  5  and  numerically.  At  first  this  seems  to  be  a 

V  M 

rather  bad  idea  since  it  should  be  more  expensive  to  double  the 
number  of  integrals.  However  this  is  not  the  case.  Since  IR  > 
»  (2ti)  exp{-lk  -R  }  both  integrals  in  (III. 7)  are  Fourier 

V  T1 

transforms  and  the  use  of  a  fast  Fourier  transform  (FFT)  algorithm 

requires  only  (N  In  N)°  operations  for  each  Integral.  Thus,  this 

leads  to  a  much  faster  algorithm  than  the  evaluation  of  the  path 

integral  by  the  matrix  multiplication  procedure.  The  additional 

efficiency  comes  from  both  the  use  of  the  split  propagator  formula^ 

(III. 5)  which  reduces  the  number  of  time  steps,  and  the  use  of  FFT 

which  reduces  the  number  of  operations  per  time  step.  The  method 

can  be  applied  equally  well  for  real  time  or  imaginary  time 
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problems.  We  emphasize  however  that  the  rapid  growth  of  labor 
with  dimensionality  confines  this  method  to  a  small  number  (S4)  of 
quantum  degrees  of  freedom;  for  imaginary  time  calculations  on 
systems  with  higher  dimensionality  Monte  Carlo  methods  should  be 
vastly  superior;  for  real  time  problems  the  Monte  Carlo  procedure 
still  has  severe  difficulties  which  take  it  out  of  contention. 

The  rela*' ionship  between  FMF  and  the  coupled  channel  method 
(CC)  depends  on  the  problem  considered  and  should  be  examined  with 
some  care,  since  a  kind  of  complementarity  exists  between  them. 
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For  strongly  quantum  degrees  of  freedom  (i.e.  those  whose 
excitation  energy  is  of  the  order  of  the  collision  energy)  the  CC 
method  is  very  efficient  since  it  requires  a  small  basis  set.  For 
weakly  quantum  variables,  however,  CC  is  very  inefficient.  The 
FMF  on  the  other  hand  is  less  sensitive  to  the  number  of  open 
states  for  each  degree  of  freedom,  and  depends  mostly  on  the 
extent  of  tlr^^e^  localization.  Localized  states  interacting 
through  localized  potentials  requires  small  grids  for  discretizing 
exp{-iTV(^) }<?|  40 ;  in  such  cases  N  is  small  and  FMF  is  very 
efficient.  At  this  point  it  should  be  apparent  why  the  use  of 
Gaussian  wave  packets  to  describe  translational  motion  would 
combine  very  well  with  a  FMF  propagation  scheme:  (a)  the  initial 
Gaussian  packet  is  spatially  localized  and  this  makes  FMF 
efficient;  (b)  for  reasonably  brief  collisions  the  packet  may 
evolve  into  a  non-Gauasian  wave  function  but  it  is  likely  to  stay 
spatially  localized;  (c)  the  use  of  ^  GWP  initial  state  is  not  a 
limitation  since  the  analysis  of  the  resulting  asymptotic  state 
easily  yields  the  S  matrix  between  many  incoming  and  outgoing 
plane  waves;  (d)  if  an  incident  wave  function  i|i(x)  is  spatially 
extended  we  can  break  it  up  in  a  sum  of  pieces  {i.e.  we  can  write 
(|i(x;0)  «  Z  <^^(x;0)  where  can  be  GWP's,  for  example)  each  having 
a  smaller  support  and  thus  requiring  a  smaller  grid.  Since  FMF  is 
exact  each  piece  can  be  propagated  independently  and  the  scattered 
wave  function  can  be  exactly  rebuilt  as  the  coherent  sum  of  the 
scattered  pieces.  In  principle  this  procedure  is  ideally  suited 
for  parallel  computing. 

For  internal  degrees  of  f reedom-which  are  localized  by 

def inltion-the  FMF  is  almost  always  convenient.  Our  calculations 

of  the  evolution  of  the  Morse  ground  state  driven  by  a  laser  show 

90 

it  to  be  extremely  efficient. 

Another  great  advantage  of  the  FMF  method  is  that  it  only 
requires  the  values  of  the  potential  at  grid  points.  We  do  not 
need  to  compute  matrix  elements  between  the  potential  and  a  basis 
set  as  in  CC,  and  the  complexities  associated  with  chosing 
potential  forms  and  basis  sets  which  are  compatible  (i.e.  lead  to 
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integrals  that  are  easy  to  compute)  does  not  appear.  This  is 
particularly  important  when  the  quantum  system  interacts  with  a 
classical  many  body  system  and  the  potential  is  generated 
nximerically  by  a  Monte  Carlo  or  molecular  dynamics  procedure. 

Ill .  The  application  of  FMF  to  the  MGWP  theory  of  Scattering 
As  in  the  case  of  the  GWP  approach  to  this  problem  it  is 
convenient  to  use  the  hybrid  coupled  channel  -  GWP  approach 
embodied  in  the  wave  function  (1.2).  Since  FMF  is  exact  the 
decoupling  of  G  .,  G  .  for  a  ^  is  no  longer  an  approximation; 

OCl  p  j 

furthermore  the  wave  functions  G^^  start  by  being  Gaussians  but 

are  allowed  to  take  any  form  Imposed  by  the  dynamics.  Thus  all 

the  approximations  made  by  MGWP  are  removed. 

Since  we  use  the  wave  function  (1.2)  we  must  solve  the  wave 

equation  (II. 5).  For  each  a  the  wave  function  is  an  m~dimensional 

vector  ^  (^;t)  =  (G  .(l?;t),  ....  G  (^;t)),  where  G  .  is  Gaussian 
at  otJ  otm  /  otj. 

at  time  t®0  only.  The  potential  energy  V^j(^)  is  a  mxm  matrix. 

To  compute  the  quantity 

exp(iTV(?  ))  .  ^  ;t)  (III. 8) 

h  Ot  11 

required  in  (III. 7)  we  must  diagonalize  the  matrix  V(^^)  for  every 

grid  point  ^  .  This  must  be  done  only  once  at  the  beginning  of 

the  Iteration  scheme.  If  we  denote  by  h(K^)  the  diagonal  matrix 

having  the  elements  A(^  ),-  =  exp(-iTX.($  where  X.(^  )  is 

^  r|ij  inlj  111 

an  eigenvalue  of  V(R  ) ,  then  we  can  write 

h 

exp(-lT  (III. 9) 

n  h  h  h 

Applying  the  matrix  B  to  gives 

f^(^)  =  •  (III. 10) 


Equation  (III. 7)  requires  us  to  Fourier  transform  each  f^(^) 
separately.  When  this  is  done  we  obtain 


(III. 11) 


F,  (^  )  =  S  <ic  I  ^  >  f ,  (I  )  . 

i  v  $  v'  T|  i  n 

n 

The  multiplication  of  4>  with  A  requires  steps  per  grid  point  ^ 
and  N  m^  steps  for  the  whole  grid.  The  FFT  leading  to  (III.  11) 
requires  (N  In  N)‘^ffl  operations. 

To  complete  the  time  propagation,  for  one  time  step,  we  must 
perform 


V 


and  this  requires  (N  In  N)^m  operations.  The  total  is  w'^m^  +  2M 
(NlnN)*^  per  time  step.  Of  course  we  must  multiply  this  with  the 
number  of  'time  steps  n  and  the  number  p  of  packets  G  required  to 
construct  the  original  wave  function;  the  total  number  of 
operations  is  Pn{N‘^m^+2m(NlnN)'^} ;  this  formula  indicates  that  this 
calculation  is  feasible  on  a  fast  computer. 

It  is  interesting  to  note  the  possibility  that  a  calculation 
in  which  we  take  too  few  time  steps  and  *  too  coarse  a  spatial 
grid  might  have  some  value  since  these  approximations  cut  off  the 
fast  excitations  (i.e.  those  transitions  having  high  frequencies) 
and  the  high  momentum  components  of  the  wave  function.  It  is 
conceivable  that  this  happens  without  strong  distorsion  of  the 

lower  and  mid  frequency  and  momentum  part  of  the  wave  function. 
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Our  experimentations  with  simple  models  confirms  this 
assumption,  but  of  course  it  does  not  imply  that  this  must  a 
general  property. 

We  note  that  it  might  appear  that  a  scattering  calculation 
require  a  large  grid  since  the  FFT  subroutine  is  such  that  the 
wave  function  is  reflected  by  the  grid  boundary.  Thus  it  appears 
that  we  must  place  the  grid  edge  very  far  from  the  scattering 
center,  to  permit  the  wave  function  to  escape  completely  from  the 
Interaction  region  without  reaching  the  border  of  the  grid.  One 
can  avoid  this  by  removing  pieces  of  the  wave  function  as  soon  as 


they  emerge  from  the  interaction  region  and  before  they  reach  the 
90 

grid.  Then  the  total  scattered  wave  function  can  be 
reconstructed  from  these  pieces  with  very  little  effort.®® 
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FIGURE  CAPTIONS 


The  coordinate  system  used  in  this  paper  XYZO  is 
fixed,  with  the  XOY  plane  on  the  surface  and  the  OZ 
auces  pointing  towards  the  vacuum,  and  describes  the 
position  of  the  center  of  mass  of  the  diatomic  AB. 
The  axis  of  the  system  XYZ  are  parallel  to  those  of 
XYZ,  and  the  center  of  the  coordinate  system  is 
moving  with  the  center  of  mass  of  the  molecule,  e 
and  ^  describe  the  onculation  of  the  molecule  with 
respect  to  the  surface.  The  Interatomic  distance  is 
frozen. 
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APPENDIX 


To  derive  equations  of  motion  for  the  parameters 
and  appearing  in  each  Gaussian  wave  packet  we  start  from 

the  expression  for  the  error  E  given  in  the  text  above  the 
Equation  (ZI.7a).  We  introduce  in  the  formula  for  E  the  Eq. 

(II. 6)  for  G  .  and  take  the  time  derivatives  indicated  in  E.  This 
leads  to  a  bilinear  expression  in  and  • 

Minimizing  E  with  respect  to  the  above  variables  (i.e.  ^ »  ^oti ' 
etc.)  leads  to 


where  n  «  0,1,  and  2.  The  resulting  three  equations  can  be 
written  as 


(A.l) 


(A. 2) 


;d5{{^-S^j{t))2[{^-S^^{t))-('^j(t)  +  I  '^j{t)-*A^j(t)) 


(A. 3) 


(^-^„j(t))]  .  |G^j(^;t)|2>+ 


2 


where 


-  -  r  ''I  *  ’’.i''* 


''nlm’aj  '  ( t ) ) "(y-y^^ ( t ) ) ^ ( Z-z^^ ( t ) ) ”| (5; t ) | ^ , 


-e  .)t/fi 


(A. 4) 


;d^(5-^^j(t)2|G^j(i5;t)  1^  . 


(A. 5) 


The  equations  (A. 1-3)  represent  13  complex  equations,  although 

(e) 

(A. 3)  is  symmetric.  F  .  is  a  sum  of  matrix  elements  coupling  the 

®'J  ritrt  ^"*0^ 

State  J  to  other  rotational  states.  P  ^  and  P  .  are  the  first  and 

aj 

second  moments  of  these  elements.  M  * is  a  3x3  real  matrix  of  all 
2nd  order  moments  (l20o’otj'  ^^llO^aj 

and  imaginary  parts  of  equation  (A. 2)  yield  equations  (II. 7a)  and 
(II. 7b)  for  P  .(t)  and  R  .(t)  respectively. 

aj  aj  ^ 

The  equations  for  the  elements  of  A^j(t)  are  not  so  easily 
separable.  We  can  combine  (A.l)  with  (A. 3)  however  to  write 


M  '  ' 

/d^(^-^  {t))2  -  - ]  |G  (K;ti 

^  “J  <^000’aj  “J 


[{^-i„j(t)).  D^j(t).(5-5„j(t);l 


(A. 4) 

Hrj-(2)  (0)  ^ 

aj  ^aj  ^000  a  j 


where 


D-..(t)  =  A..^(t)  A.,(t)-  A..(t) 


3 


This  can  not  be  reduced  to  a  simple  matrix  equation  for  ( t ) . 
However,  the  matrix  of  equations  (A. 4)  is  symmetric,  and  three  of 
the  equations  are  redundant.  Thus,  we  solve  for  a  vector 

containing  the  six  non-redundant  elements; 

Upon  Integration,  and  a  bit  of  matrix  alegbra,  equation  (A. 4)  can 
be  written  as 


-t  +  c  .  =  0, 
aj  aj  aj 


where 


,(0) 


■J(2) 

^^OOO^aJ 


T  .  is  a  vector  of  all  second  order  moments  of  G  j{^:t); 


Thus 


and  we  arrive  at  matrix  equation  (II. 7c),  where 


Finally,  using  the  fact  that 


(^-^^j(t)) 


oj  aj' 


as  equation  {II. 7d)  for  y  j(t). 


we  can  write  equation  (A.l 

As  mentioned  in  the  text,  the  ability  to  compute  these  three 
dimensional  moment  and  potential  Integrals  analytically  is  a 
tremendous  advantage  of  using  the  Gaussian  basis.  The  moments  are 
straightforward  to  compute  and  can  be  derived  from 


<  ^OOO^oJ 


/d5 


-  I  tIm(*AtS^^]  -|lm(y^,(t)) 


otj 


aj 


8det(Im(  A^j(t))) 


All  moments  can  be  generated  by  talcing  derivatives  of  (!„««)  4 

UUU  Ocj 

with  respect  to  the  elements  of  A^j{t).  For  example 


^^200*aj  “ 


h 

2 


3 (  A  .  ( t)  )  ) 

a  j  XX 


‘■000 


^^OOO^otj 


"”<^ai<")vvI">(A«i(t))„-(Im(A^.(t))^v 

detdml’^^it)  )  ) 


All  odd  order  moments  are  zero. 

The  potential  integrals  in  F^j  can  all  be  written  as  sums  of 
terms  of  the  form 


c  ;  d5  (5-i?^j(t))"  e~^’  ^  , 

where  c  is  some  collection  of  constants,  and  the  wave  packet  and 
potential  parameters  are  contained  in  A  and  ?.  For  n  =  0,  we  can 


